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An overall rise in fuel prices has led to an increasing 
interest in the design of autopilots for ships. The purpose 
of the automatic steering control is to minimize the prepul- 
sion losses, which are caused by added drag due to steering 
of the ship. Minimizing a performance criterion based on 
added drag due to steering can reduce fuel consumption. 
Claims by many researchers indicate that a carefully 
designed ccntroller could save from one to two percent of 
fuel. For large containershirps this could amount to more 
than $100,000.00 per year savings. 

To study the optimization problem, models of toth the 
Ship and its operating environment are reguired. what type 
of ccmputer model should be used to represent the ship? 
Chapter two addresses the development of several nodels. 
Since the best model was desired it was decided to use the 
equations of motion to simulate the Ship in our Fortran 
progran. The basic Nomoto models give an adeguate descrip- 
tion of ship steering dynamics for design. The Ncmoto 
second- and third-order models were developed from the egqua- 
tions of motion as defined by a series expansion including 
all terms (both linear and nonlinear) for which nydrodynapic 
coefficients were available. An interactive program that 
utilized the Nomoto mcdels to model the ship was also used. 
Two independent programs were developed to aid in the design 
of the ccntroller. 

What 1S an adequate cost function which represents the 
added drag due to steering? Chapter three addresses’ the 
classical cost functicn used by many researchers. 

Since a variety cf control algorithms are possible one 


must ask if one algorithm provides a lower minimal cost than 
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another. Chapters fceur, five and six address the selection 
of the ccntroller which provides the minimum value of added 
drag due to steering. 

Ship dynamics change with operating conditions such as 
ship speed , sea state , and encounter angle. Therefore an 
adaptive controller must be used to provide minimum added 
drag due to steering. Chapter seven development of an 
approach to an adafgtive controller utilizing satellite 
infornaticn. 

Conclusions were drawn from these experiments and are 
presented in Chapter eight. This thesis investigated only 
course keeping with emphasis on minimizing rudder and yawing 
activity to reduce fuel consumption. Presented in this 
Chapter are recommendations for future study where the 
objective is track fcllowing which would be important for 
ships reguired to follow stringent routes. It 1s alsc impor- 
tant fcr other systems such as. satellites, nissiles, 
aircraft, where the controller minimizes yaw error to keep 
the system on track. 
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II. COMPUTER MODELS 


The model which best represents ship-Ssteering dynamics 
is a Taylor's series expansion of the force and moment rela- 
tionships around a selected steady-state operating point. 
The resulting equaticns are commonly known as the eguations 
of motion [Ref. 1]... A computer program was developed using 
known available data on the hydrodynamic coefficients for 
the SI-7 containership to provide a computer simulation of 
the ship. The computer program is shown in Appendix A. 
Figure 2.1 shows the block diagran. Small yaw command 
angles are used, for example YAWC= 1.0 / 57.296 represents a 


yaw ccaumand change of one degree. 






EQUATIONS 
OF 
MOTION 





CONTROLLER 





ee i ee Sa — 
-—- - ee te 


Pigure 2.1 BLOCK DIAGRAM 
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To oktain the Noszoto second- and third-order transfer 
functions from the equations of motion, the function mini- 
mization subroutine was used to obtain the coefficients. 
Figure 2.2 shows the scheme used to obtain the Ncsoto 
transfer functions. The computer program is shown in 


Appendix A. 











FUNCTION 
MINIMIZATION 
SUBROUTINE 


ae i 2 
Jef yw dt 





Figure 2.2 DETERAINATION OF NOMOTO MODELS 


The Nomoto models were checked against analytic results 
from linearized equations. 
Proceeding to the second-order Nomoto equation: 
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¥(S)/5(S) = K/S* (1+1*S) (2.1) 


Deriving the second-crder Nomoto transfer function from the 
yaw equation only, the result is 

¥{S)/6(S) = 0.040893/S* (1+8.539932*S) 
and using function mirimization as in Figure 2.2 

V¥(S)/6(S) = 0.0409221/S* (1+8.55207382*S) 
and the agreement is obvious. Using function mininization 
with both yaw and sway equations with linear terms only, the 
results are: 

¥(S)76(S) = 0.1072741/S* (143 1.9199524* S) 
If the nenlinear terms are included but the perturbaticn is 
snall 

W(S)78(S) = 0.1072082/S* (1431.8907013*S) 
and it is clear that the nonlinear terms contribute little. 


Proceeding to the third-order Nomoto equation: 


V(S)/5(S) = K*¥(1#T24S) /S* (14#7P1¥*S) * (14TP2*S) (Oe) | 


The parameters were calculated and checked by using function 
Minimization as in Figure 2.2. The results are given in 
Table 1. It is clear that the answers obtained by function 


minimization agree clcsely with the analytic solutions. 


TABLE 1 
THIRD-ORDER NOMOTO MODEL FOR THE SL-7 


speed K TZ TP1 TP2 

knots calc comp calc comp calc comp calc comp 
160 «6.0738 .0738 22.57 22.95 12.986 12.946 107.583 107.583 
23 21067 1.1061 15.67 15.70 9.014 9.006 75.130 74.846 
32 1477 1477 11.28 11.28 6.470 6.467 53.793 53.793 
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Analytical eguations used to calculate seccnd-crder 


Nomotc transfer functicn coefficients are: 


K =N N ° T =-({(IT -N. N 
6 4 Yr F (T r ys ic 


Analytical eguations used to calculate third-crder 


transfer function coefficients are: 


K = (No *Y,/¥)/ (NN *(Y_-40) /Y_) 

TZ =-((U-¥,) #Ny -N-*Y,)/(Y #0, -N FY. ) 

TPIFIP2=— ( (M-¥2)* (TN 2 )-NS*Y 2) / (WN (Y -M4U)-Y #8) 

TP 1+TP2=((N-Y-) #N + (1,-N2) *Y +N, *(¥ -4#U) +¥2*N ) 
/(NS*(Y NSO} Yea} 


The nondimensionalized hydrodynamic coefficients for the 


SL-7 ccntainership are shown in Table 2. 


\ TABLE 2 
HYDRODYNABIC COEFFICIENTS FOR THE SL-7 


axial force lateral force moment z-axis 
x" = -0.0001 y? = -0.00758 N° = -0.00213 
V V 
x? = -0.9003 y? = 9.0023 Ne = -0.00105 
uu r i 
eee = 90.0039 ¥". = 9.00145 Ne = -0.0007 
x! = -0.0012 y? = 0204 N! = —-0.015 
VV Waige Vie 
x? = -0.0005 y? = -0.008 ne = -0.008 
66 vrr V Ge 
y! = -0.03 N? = 9.01 
VVV VVV 
Ma 0.003 N' p= 70-006 
Yess = -0.0005 Ss sun 0.0001 
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eS Es =P ea eS 


In recent years, many have studied the protlem of 
{[Ref. 2] [Ref. 3] [Ref. 4] [Ref. 5] [Ref. G6] [Ref. 7} 
[Ref. 8] [Ref. 9] [Ref. 10] [Ref 11] [Ref- 12] [Ref. 13] 
{Ref. 14] optimizing an automatic ship-steering contrcller 
for minintun fuel consumption. It is well known that addi- 
tional drag is intreduced by steering and that both the 


rudder moticn and the yawing motion contribute to this added 


drag. A measure of the added drag given as a cost function 
is 
ae 
me Aste ( At. 2 + 62) dt (sel) 
2 


where W, = yaw error 


or 
i 


rudder angle 


>” 
il 


weichting factor 

While this expression 1S an approximation, it is conven- 
jent for shipboard use because y, and 6 are readily measur- 
able. There is no general agreement on numerical values for 
the weighting factor, X\ , and in this study the valves used 
were chosed from the work of R.E. Reid [Ref. 7] <or the 
SL~7. 

The weighting factors for the operating range of the 
ship are shown in Tatle 3. 
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TABLE 3 
WEIGHTING FACTOR 


ship speed weighting factor 
(knots) 
16 16.796 
23 8.128 
ae 4.2 


Reid's work shows the relationship of weighting factor 
to the closed-loop natural freguency, mass, pivetb pein, 


ship speed, xX" and X'.. hydrodynamic coefficients. It is 


a 
shown in cunts ae SG Reid chose a closed-loop natural 
freguency of 0.05 rad/sec which experimentally shcewed at 
this frequency, the weighting factor in the cost function, 
provided good representation of the added drag due to 


steering. 


N= 24M#(14X" ) # (CP/L) #02 /(P/2)* (14K. #02) (3.2) 
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IV. CONTHOLLER DESIGN USING ICSOS 


The Interactive Control Systen Optimization and 
Simulation (iCSOS) rackage finds optimum values for unknown 
{free) parameters in a control system design problen and/or 
performs situlation of the system. An example of usage of 
ICSOS is shown in Appendix B. 

In preliminary studies ICSOS was used with Nomotc models 
to study ccntroller characteristics in calm water. The func- 
tion Binimization sukroutine adjusted the controller paranme- 
ters to minimize the cost function. Figure 4.1 shows the 


scheme used to evaluate the controller. parameters. 








NOMOTO 
MODEL 






Porometers 


FUNCTION 
MINIMIZATION 
SUBROUTINE 


1 
Jef (+8041 
0 e 












Figure 4.1 OPTIMIZATION OF CONTROLLER 
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Reid [{Ref. 7] uses the second-order Nomoto model of 
equation 2.1 for the Si-7 and also uses a cCcntrcller 


descrited by 


Go(S) = K1*(1#T1*S) ~ (1#T2*S) (4.1) 


His results are given in Table 4. 


TABLE 4 
REID'S RESULTS 


speed plant weighting controller gains 

knots K a factor K 1 T1 je 
16 0.1084 90.36 16.796 064555 Bea.) | VO sore 
23 6.1556 64.67 8.128 0.37609" “62260 8.308 
32 Ve 2e7 e854 5 : 0.3188 44.92 7.066 


Using this plant and weighting factor values but applying 


IcSOS, results were chtained and Shown on Table 5. 


TABLE 5 
ICSOS RESULTS 


speed plant weighting controller gains cost 

knots K ak factor K1 T1 2 J min 
16 - 1084 90.36 16.796 .454616 90.3459 10.0215 340.864 
23 -1556 64.67 8.128 .373171 64.0658 8.4640 139.9916 
32 02167 45.45 : 318645 45.4475 7.0662 60.828 


In each case the controller zero (1/T1) cancels the 
plant pele (1/1). Additional experiments consisting of 
inserting arbitrary numbers in the Nomoto equation and 
repeating the computer run indicated that this will always 
be true. That is, to minimize the cost the plant pole is 
cancelled and a new fole location determined with apfpropri- 


ately adjusted gain. 
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The simple contrcller of Eguation 4.1 41s an arbitrarily 
chosen structure. To determine the effects of more complex 
controllers three additional structures were chosen as shown 
_ain Figure 4.2. Each of these was used with the Ncmoto 
second-order model for the ship at each of the indicated 
speeds. The results are shown in Tables 6, 7, and 8. 
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Figure 4.2 VARIOUS STRUCTURES FOR CONTROLLERS 


These results are very interesting. At 16 knets’ the 
contrcller gain {K1), controller zero (1/T1) and ccntroller 
pole {1/12) are essentially the same for all structures. For 
structure B, which includes an additional pole, the function 
Minimization subroutine tries to drive the additional pole 
to infinity, and no doubt would have done so if the calcula- 
tions had continued. Por structure C, which has two foles 
and two zeros, a zero and pole cancel indicating that they 
are not needed or wanted. For structure D, the integrator 
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TABLE 6 
SIBULATION RESULTS = STCADY SPA rt eeolvc yo ce. 


CALM WATER FOR VARIOUS CONTROLLERS 
FOR FIXEE SHIP SPEED { 16 KNO s_} 
BOMCTO SECOND-ORDER MODEL (K=.1084 tr. 0. 36) 
Y= 16- OPTIMAL PARAMETER GAINS O 
VARIOUS CCNTROLLERS , COST FUNCTION. 
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TABLE 7 
SIMULATION RESULTS —- STEADY STATE 600 SECS 


CALM WATER FOR VARIOUS Se Ca 
FOR FIXED SHIP SPEED 23 KNOTS ) 
NOMOTO SECONT-ORDER MODEL ({K=. 1556,T=64.6/7) 
A= 8. OPTIMAL PARAMETER ChINS OF 
VARIOUS CCNTROLLERS » COST FUNCTION 
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TABLE 8 
SIMULATION RESULTS — STEADY STATE 600 SECS 
CALM WATER FOR VARIOUS CONTROLLERS 
POR FIXED SHIP SPEED ( 32 KNOTS s_) 
NOMOTO SECOND-ORDER MODEL (K=.2167,T=45. 45) 
OPTIMAL PARAMETER 6kINS OF 
VARIOUS CONTROLLERS , COST FUNCTION 
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gain 1S driven to zero. The same pattern of results is 
obtained at 23 knots and 32 knots. Note that in all cases 
the minimum cost is essentially the sane, as would be 
expected since all ccntrollers are the same. 

Using the computer method of Figure 4.1 and the Nemoto 
third-order models of Table 1, controllers A, B, C of Figure 
W.2 were optimized... The results are shown in Tables 9, 103 
and 11. 


TABLE 9 
SIMULATION RESULTS - STEADY STATE 600 SECS 
CALM WATER FOR VARIOUS CONTROLLERS 
FOR FIXED SHIP SPEED aa oo ) 
OMOTO R 
(K=. oe” 2a TZ=22.56/3,TP1=1 D8 458, TP 2= 107.5853) 


16.796 , OPTIMAL PARAMETER GAINS OF 
 TARTOS CONTROLLERS , COST FUNCTION 


contr ccentroller gains cost 
K 1 T 1 r3 T4 J nin 
A Q0.648486104 90.C€S94 15.27712 - - 3704023 
B 0.6441367 84.826 15.78691 .24598 os 374.3808 
C 026151139 107.5782 8.73520 12.9368 24.9676 369.9297 
TABLE 10 
SIMULATION RESULTS - STEADY STATE 600 SECS 
CALM WATER FOR VARIOUS CONTROLLERS 
FOR FIXED SHIP SPEED 23 KNOTS ) 
HO ACTO Cee de R MODEL 
iK4 1067, TZ=15. 675,TP1=9.014, TP2=75. 13 
8.128 , OPTIMAL PARAMETER GAINS 0 
VARIOUS CENTROLLERS , COST FUNCTION 
contr ccntroller gains cost 
K1 T1 T T3 T4 J min 
he O.5224258 63.13609 12./2212 = - 152.2920 
B QO.5216467 64.93709 12.63218 .0505174 = fe 2.5303 
C (€.5001907 75.14&52 6.527490 9.039928 18.260 152.2800 


Of zajcr interest 1s the fact that the difference in 


"cost" between A, B, C is less than one per cent. At each 
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TABLE 11 
SIMULATION RESULTS —- STEADY STATE 600 SECS 


CALM WATER FOR VARIOUS CONTROLLERS 
FOR FIXED SHIP SPEED { 32 KNOTS ) 
CTO THI RD-ORD R MOD 
(K=. 14771 rae 11. 2833 TP1=6. 4699, 132-53. = 7931) 
Y= MAL PARSER ETER GAI 
VARIOUS CONTROLLERS , COST FUHCTION 


contr contrcller gains cost 
K1 Tl 2 nS T4 J min 

A 0.427633 48.66C48 10. 744 8&5 - - 68.09039 

B 0.298732 89.40696 15.01033 .0597786 - 69.3 72355 

C 0.417991 53.6961 4.970016 6.294354 13.85724 68.04735 


speed (16€,23,32 knots) controller Cis "BEST", but the 
difference is slight. Examining the parameter values 
obtained fcr contrceclliec C, 1t is seen that at all three 
speeds beth poles of the ship are essentially cancelled by 
zeros of the controller. 

These results seem to indicate that the dynamics cf the 
plant determines the cptimum structure for the contrcller. 

Using a state-feedback controller and Nomoto third-crder 
wodels of Table 1, the contrceller was optimized for varicus 
ship speeds. Figure 4.3 shows the scheme used to evaluate 
the state-feedback ccntrollec. 

Using the scheme of Figure 2.2, with no change in cost 
functicn or weighting, the optimal gains and costs’ were 
deterfined as shown in Table 12. 

When ccmparing the state-feedback controller with 
contrcller C, it is seen that at each speed controller C has 
a lower cost. Among the controllers consired, controller C 
is "BEST" when using the Nomoto third-order model, although 


the differences in ccst are not dramatic. 
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Figure 4.3 OPTISIZATION OF STATE PEEDBACK CONTROLLER 


TABLE 12 
SIBULATION RESULTS — STEADY STATE 600 SECS 
CALA WATER FOR VARIOUS SHIP SPEEDS, 


OPTIMAL PARAMETER GAINS FOR 
STATE~FEEDBACK CONTROLLER 


speed Nemoto third-crder plant weighting ccntroller ccst 

knots K TZ TP TP2 factor K1 K2 J min 
16 .0738 22.567 12.946 107.583 16.796 4.426 78.004 370.711 
23 ©1067 15.675 9.01% 75.13 2128 3.103 45.649 152.596 
32 .1477 11.283 6.470 53.793 oa 0240 27.896 68.2513 
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Vo CONTROLLER DESIGN USING FORTRAN PROGRAMS 





The Fortran program referenced in Chapter two which 
provided a computer Simulation of the SL-7 ship was rodi- 
fied. A function minimization subroutine was coulped to the 
Simulaticn and used the subroutine to adjust ccntrcller 
parameters to minimize the cost function and to evaluate the 
Binimum cost. Figure 5.1 shows the scheme used to evaluate 
the contrcller parameters. This program was used for conpar- 
ison to ICSOS. The computer program is shown in Appendix A. 
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Figure 5.1 OPTIMIZATION OF CONTROLLER USING FORTRAN PROGRAM 
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Using the computer method of Figure 5.1 and the nenli- 
Near €guations of motion, controllers A, B, C of Figure 4.2 
were optimized. The results are shown in Tables 13, 14, and 
a. 


TABLE 13 
SIMULATION RESULTS - STEADY STATE 600 SECS 


CALM WATER FOR VARIOUS CONTROLLERS 
FOR FIXED SHIP SPEED § 16 KNOTS ) 
ECUATIONS OF MOTION 
A= 16.796 , OPTIMAL PARAMETER GAINS OF 
VARTOUS CCNTROLLERS , COST FUNCTION 


contr contrecller gains cost 
K 1 lt 2 is T4 J mln 
A ~-€4€401 89.81704 15.381699 = - 1.128189 
B sozeus0 9C.67294 159542297 0.92091336 = o t73e235 
c Same peewee 8.597099 132353928 25.21362 1.126307 


TABLE 174 
SIMULATION RESULTS —- STEADY STATE 600 SECS 


CALM WATER FOR VARIOUS CONTROLLERS 


ECUAT MOTION 
OPTIMAL PARAMETER GAINS OF 


A= 8.128 
VARIOUS CCNTROLLERS , COST FUNCTION 
contr controller gains cost 
K1 T1 T2 T 3 TY J min 
me). 522106 66233122 12.83327 ~ = 0.4640879 
B 2.865869 66.151£z 13.01183 0.92783 = 024857854 
C .503967 74.79771 6.65880 9.20533 18.4022064 0.46 36095 


These results agree with those obtained by ICSOS and 
contrcller C provides the minimum cost. 

If the assumption that the steering dynamics of the ship 
is adequately modeled as a second-order system is valid, 
then cnly two states are needed for feedback. Pot avthird— 
order system three states are required. The ccntroller 


structures are shown cn Figure 5.2 and 5.3. Using the scheme 
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TABLE 15 
SIMULATION RESULTS ~ STEADY STATE 600 SECS 


CALS WATER FOR Seo CONTROLLERS 
FOR FIXEC SHIP SPE D eee 32 i oe ) 
EQUATIONS OF M 
A= 4.2 PTIMAL PARAMETER GAINS OF 
VARIOUS CCNTROLLERS g COST FUNCTION 


contr controller gains cost 
K1 T1 T2 73 TY J min 
A 2428804 48.65540 10.874426 = = 0. 2072417 
B .298732 89.40696 15.010330 0.01 0.2118334 
C .417333 53.09654 5.096548 6.474857 14. 0205 0.2071124 


of Figure 5.1 with no change in cost function or weighting, 
the optimal gains and cost were determined as shown in 
Tables 16 and 17. Comparing costs, there is little differ- 
ence between the twe state system andthe three _ state 
systen. The comparison between state feedback with 
contrcller C, it is seen that at each speed controller C is 
Fetter, but not much Letter. 
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Figure 5.2 TWO STATE SYSTES 
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Figure 5.3 THREE STATE SYSTEM 


: TABLE 16 
SIMULATION RESULTS - STEADY STATE 600 SECS 


CALS WATER FOR VARIOUS SHIP SPEEDS 
EQUATIONS OF MOTION 
OPTIAAL PARAMETER GAINS FOR 
THO STATE SYSTES 


- speed gains weighting cost 
knots K1 K2 . factor J min 
16 4&.4033689 7725041656 16.796 1.128771 
23 3.9889006 45. 2637787 8.128 - 4646050 
; 202342062 2726808014 - 2075207 


32 
| / 

Note that for the Nomoto nodel studies yaw error and 
rudder angles were measured in degrees; when the eguations 
of mction were sinulated yaw error and rudder were in 
radians. Thus the numerical values of the cost, J, are 
different. y 

Transient response plots for controllers A, 5B, C, and 
three state-feedback at Ship speed 32 knots are shown in 
Figures 5.4 - 5.9. . 
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TABLE 17 
SIMSULATION RESULTS ~ STEADY STATE 600 SECS 


CALA WATER FOR VARIOUS SHIP SPEEDS 
EQUATIONS OF MOTION 
OPTIA PARAMETER GAINS FOR 
THREE STATE SYSTEM 


speed gains weighting cost 

knots K1 K2 K3 factor J pin 
16 4.8617249 87.7073364 99.9802704 16.796 1.128289 
Ze 326630983 56.2784882 88.5913391 8.128 4643548 
32 205967150 33.27571444 41.3186035 4.2 © 2074225 
co FN ee ~~ ane 
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Figure 5.4 YAW ws. TIHE (controller A, B, C and state-feedback) 
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Figure 5.5 RUDDER ws. TIME {controller A, B, C and state-feedback) 
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Figure 5.6 YAW AND RUDDER vs. TIME {(controlier A) 
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Figure 5.7 YAW AND RUDDEB vs. TIME (controller B) 
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Figure 5.8 YAW AND RUDDER vs. TIME (controller C) 
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Figure 5.9 YAW AND RUDDER vs. TIME (state-feedback ccntroller) 
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VI. CONTEOLLER DESIGN IN SEA STATE 

A sea state generator was coupled to the Fortran 
program, so that the function Minimization subroutine could 
be used in the presence of the sea state. The sea state 
generator was an elaborate program obtained from DINSRDC. 
This pregram generates added mass and inertia values as 
functicns of encounter freguency and also calculates the 
forces and moments. The forces and momentS are generated 
and stored in a look up table which was coupled to the equa- 
tions cf motion. Figure 6.1 shows the scheme used to eval- 
uate the controller parameters. The computer prograr is 
Shown in Appendix A. 

The cptimal gains obtained by the calm water study of 
Chapter five were use aS the initial guess in evaluating the 
optimal ccntroller parameters in the presence of a Seaway. 
For ccmparison, studies of the value of the cost function 
uSing calm water gains in sea state were obtained; then the 
function finimizaticn subroutine was allowed to adjust 
contrclier parameters in the presence of several sea states 
and encounter angles. The entire study was done at a ship 
speed of 32 knots. The added mass and inertia change with 
tespect to encounter frequency as shown in Figures 6.2 and 
6.3. Figure 6.3 is ncndimensionalized by dividing the added 
inertia by the mass of the displaced water and the sguare of 
the length ketween perpendiculars. To convert back to dimen- 
Sionalized units of lb-ft-sec2, multiply the graph pcints by 
2-581E12. Since the sea state 1s represented by irregular 
waves, the waves impinging on the ship hull contain the 
total energy density spectrum composed of many freguencies 
and the ship responds to an average value of added rass and 


inertia. The values used for this study was obtained at 
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encounter frequency cf 0.75 rad/sec from our sea _ state 
generator. This freguency gave us values for added tass and 
inertia representative of an average value. The energy 
density spectra for various sea states are shown in Figure 
6.4. The added mass for Sway was changed from 2.6457E06 
lb-sec2/ft for calm water to 2.3043E06 Il1b-seceyft for a 


seaway. The added moment of inertia for yaw was changed 
from 1.42E11 il1b-ft-secé Cor calm water to 1.5096E11 
lb-ft-sec® for a s€away. All other hydrodynamic cceffi- 


clients were kept constant at calm water values. The results 
are shown in Tables 18 - 25. In certain sea states and 
encounter angles the calm water optimal gains performed well 
aS Shown by calm water cost value when compared to sea state 
cost value. In most cases the function minimization subrou- 
tine fcund new gains with lower cost function values in 
seaway as compared tc using calm water gains. In the caln 
water evaluation, the system was perturbed with a one degree 
course chanye, but the course change was not included in the 
seaway tests. The difference in cost values is attrikuted to 
the difference in operating conditions. 

Using the Proportional, Integral and Derivative (PID) 
contrcller Equation €.1 with no change in cost. function , 
the function minimization subroutine was used to adjust 
contrciler parameters to minimize the cost function and 
evaluate the minimum cost. The results are shown in Table 
26. When comparing the PID with controller A, it is seen 
that at each enccunter angle, ccntroller A 1s better. These 
results agree that in a seaway controller A frovides the 
minimum cost. , 

Table 27 shows ccmparison of the minimum cost function 
for ccentroller A, controller C, and PID. The study was done 
at ship speed of 32 knots and at sea state 4. Controller A 
provides the minimum cost. 
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-Pigure 6.1 OPTIMIZATION OF CONTROLLER IN SEA STATE 


The optimal gains obtained in the presence of sea state 
was dcne over using a Simulating time of 600 secends. the 
sea state program is designed to provide gradual increase in 
the forces and moments during an initial time interval.. 
This is done to minimize initial condition transients in the 
Ship dynamics. There will-unavoidably be some transient 
effects, however, and these could affect the value of the 
cost, J, determined during the 600 seconds of sinulation. To 
determine whether such initial transients had. any signifi- 
cant ccentribution to the value of J, additional simulation 
runs were made with the controller parameters fixed at their 
optimal values. However, evaluation of the cost, J, was 
started only after 300 seconds of simulation had elapsed. 
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Figure 6.2 ADDED MASS vs. ENCOUNTER FREQUENCY 
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Figure 6.3 ADDED INERTIA vs. ENCOUNTER FREQUENCY 
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TABLE 26 
SIMOLATION RESULTS - STEADY STATE 600 SECS 


FIXED SHIP SPEED (32 KNOTS) IN A SEAWAY 
SHIP MCUEL: Pune OF MOTION 


SEA STATE 4 
PID CONTROLLER 


encounter centroller gains sea state 
angle KI Tl T J min 
30 95263100 4.20720660 2635606705 8029 Coo 
69 -68631890 12.5794449 8.2121658 .09730512 
90 225809155 12.4257) 589 S77 G1oe CO to 
120 4.9198265 12.5986176 .6/7592590 5919. 708960 
150 12.3970&23 15.7682953 Vos loot. 4 2 ae 
8(S)/¥, (S) = K1 + K1*T1*S / (1+#T2%*S) **2 (6295 
TABLE 27 


COMPARISON OF TEE MINIMUM COST 
SHIP SPEED (32 KNOTS) 
SEA STATE & 


enccunter controller controller controller 
angle PI 


degree J min J min J min 
30 -022854677 - 034033 ~02985619 
60 7093756 97 098914 - 09730512 
90 1.51713 40 15/206 1359 15950 
120 9.9917300 10.3530 On 08980 
150 16.67052ZC 2023956 17.427290 


The value obtained was then doubled and compared with the 
result of evaluating J over the full 600 seconds. Corparison 
of Table 28 with cost values in Tables 18, 19, 20, 21 shows 
only small differences. 

Io oktain insight into the stochastic process of irreg- 


ular seas, a deterministic process was studied. The Fortran 
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TABLE 28 
EFFECTS DUE TO TRANSIENT AND GRADUAL BUILD UP OF SEA STATE 


INTEGRATION OF COST FUNCTION ( 300 TO 600 SECS) 
FIXED SHIP SPEED {32 KNOTS) IN A SEAWAY 


sea encounter ccntroller gains cost ecst 

state angle K 1 T1 2 J min 2*J min 
1 60 12.403329& 10.650075 1.1086683 .0641122 .1282244 
2 E0 moe Do oe eco, Joloos 2. 30770530 1997S) .0399962 
4 60 -62012090 40.80556€0 19.606873 .0515974 .1031948 
6 60 127228041 8.401714740 .51471250 .7906179 1.581236 


program was modified to minimize the cost function in the 
presence of a regular sea. To allow compariscn with 
previcus work the encounter frequency of 0.75 radysec was 
used and scaled the amplitude of the regular sea to its 
Frosrective sea state. The entire study was done at a ship 


speed of 32 xnots. The results are shown in Tables 29.530 4 


Table 29 shows that for regular seas the ccntrcller 
parameters do not ckange significantly for different sea 
states; but as sea state increases, the cost value increases 
due to the increase in yaw moment and sway force on the 
Ship. Tables 30 and 31 also show that the controller parame- 
ters donot change Significantly from sea state to sea 
Staaec « However, an encounter angle of 90 degrees Shows a 
relatively high cost compared to costs calculated for 60 and 
120 degrees at a given sea state. To account for this 
anomaly, the following is suggested. In the regular sea, 
the added mass and inertia were known for a given encounter 
frequency, while in the irregular sea a representive average 
value was used. The rethod used to obtain the average Bight 
not represent the actual average. Also, it seems reascnable 
to suppose, that the assumpticns of the function weighting 
factor are satisfied for all encounter angles; that is, the 


weighting functicn (Eg. 3.2), which appears in the cost 
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function (Eg. 3.1), does totally represent the added drag 
for all encounter angles. Future study is needed tc answer 
these guestions. 

The sea state in the deterministic model is represented 
by cegular waves. On this description, the waves impinging 
on the ship hull ccrrespond to only one freguency in the 
energy density spectrum. In the case of irregular seas, 
however, the spectral components change for different 
states, aS shown in Figure 6.4. Thus comparison of the 
contrclléer parameters obtained for regular seas with results 
for irregular seas is not justified. The function tinitiza- 
tion subroutine adjusted contrcller parameters to minimize 
the cost function fcr either case (irregular or regular 


seas) as shown in tatles 32 and 33. 


TABLE 29 
SIMULATION RESULTS - STEADY STATE 600 S#CS 
FIXED SHIP SPEED (32 KNOTS) IN A REGULAR SEAWAY 
SHIP MCDEL: ee OF MOTION 


C=0.0 
ENCOUNTER FREQUENCY = 0.75 RAD/SEC 
ENCOUNTER ANGLE = 60 DEGREES 


CONTROLLER A 


sea ccentroller gains cost 
State K 1 Ta lz J min 
1 - 1449795 141. 383179 32.9405670 -000764582 
Zz 21534657 129 20 7a 31.4042358 003056434 
4 - 1514665 1359-9 e7 ay 32.9749756 -009345479 
6 - 1533340 135.488495 33.2869 532 022174600 


Note that in Ecth the deterministic and stochastic 
models, among the ccntrollers considered, controller A is 
"BEST" ian a seaway disturbance, although the differences in 
cost are not dramatic. 

Finally, the observed dependence of optimal ccntroller 
gains on sea state and encounter angle suggests that an 
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TABLE 30 
SrIPULATITON RESULTS — STEADY STATE 600 SECS 
FIXED SHIP SPEED (32 KNOTS) IN A REGULAR SEAWAY 
SHIP MCDEL: ae OF MOTION 
» ENCOUNTER Bee aaa a = Pi RAD/SEC 


STATE 
CONTROLLER A 
encounter ccntroller gains cost 
angle K 1 iL 2 J min 
30 223517043 1022021975 Ze as Ji oer me 2595300 
60 - 1514665 liga. 7987 38 S27 0 BUO9SHSETS 
90 - 4968442 66.546493 49.7598267 2048143090 
120 Soe 250 149.540543 33.6013489 038937880 
150 4536914 Fee 506523 Saco oo aro 2534 153 
TABLE 31 
SIMULATION RESULTS - STEADY STATE 600 SECS 
FIXED SHIP SPEED (32 KNOTS) IN A REGULAR SEAWAY 
SHIP MCDEL: Ween OF MOTION 
ENCOUNTER FREQUENCY = 0.75 RAD/SEC 
SEA STATE 6 
CONTROLLER A 
encounter ccntroller gains cost 
angle K1 T1 T2 J min 
30 mee TOO Z 100.122940 28.0581207 mo 709767 o 
6G - 1533340 135.488495 a oOo Cc 29022174600 
90 =D 21050 7 S2anl oD 30 49.9858093 wile? 72580 
120 - 1414837 142.695160 85.317 ise -091541€50 
150 4587426 71.451385 33.4558024 - 148615825 


adaptive centroller must be used to provide a continuous 
Minigum on the cost function. 

After obtaining tke optimal gains for controller A, to 
observe the behavior of the rudder and yaw motion of the 
Ship, transient resfcnse plots were obtained for controller 
A at ship speed of 32 knots and sea state 4 for various 


encounter angles as shown in Figures 6.5 - 6.14. Note the 
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REGULAR SEAS CONTROLLER GAINS 
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as the encounter 


rudder and yaw amplitude 


angle increased. 


the increase in yaw moment 


This is due to 


and sway force on the ship. 
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Figure 6.7 YAW vs. TIME 60 DEGREES 
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VII. AN ADAPTIVE CONTROLLER 

In a s€away, the contrcller gains changed dramatically 
for changes in sea state and encounter angle. An adaftive 
contrcller must be used to provide continuous operation on a 
minimum cf the cost function. This Chapter addresses a 
theoretical design cf an adaptive controller. 

In the future, there will ke better measurement cf navi- 
gation than can be fgrovided by conventional equipment on 
koard a ship. Presently the Navy 1S involved in a program 
that will provide precision Navigation data. The 
NAVSTAR/GICBAL POSITION SYSTEM (GPS) [Ref. 15] [Ref. 16] 
[Ref. 17] will provide extremely accurate three-dimensional 
positicn and velocity information to users anywhere in the 
world. The position determinations are based on the measure- 
ment of the transit time of RF signals from four satellites 
of a total constellation of eighteen. This system is sched- 
uled to be fully oférational in 1988. At present (1984) 
there are four NAVSTAR/GPS satellites in operation which 
allows three to fcur hours per day of navigaticn .tize. 
Already the Texas Instrument Company markets a receiver for 
this System where GPS can be used. } 

The Navy Remote Ccean Sensing System (NROSS) [Ref. 18] 
will te able to determine wind velocities over the world's 
oceans with an accuracy sufficient to determine ccean 
surface waves. It's objective will be to acguire glotal 
ocean data for operation and research use by both the #ili- 
tary and civil sectors. This system is scheduled to launch 
its first satellite in June 1989. 

The scheme fcr an adaptive controller is shown in Figure 
en |= Having stored the optimal controller parameters ina 
look up table as functions of ship Speed, sea state, and 
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encounter angle, the ship operating condition must Le known 
so that the table is useful. NAVSTAR/GPS would identify 
ship speed and NROSS would identify sea state and encounter 
angle. The optimal farameters can then be looked uf and 
inserted into the controller. This should place system cper- 
ation fear the minimus J. To ensure fine tuning, a micro- 
processor programmed, with the function minimization or-line 
in machine language, with inputs of yaw error and rudder 
moticn of the ship would accomplish the fine tuning rapidly. 
Since tke subroutine is written in Fortran (as used for this 
study) this would be inappropriate for on-line use. 

The adaptive controller can be performed with digital 
circuits rather than analog components. Garcia [Ref. 19] 
demonstrates the process for converting an analog ccntrciler 
into a digital controller. Figure 7.2 illustrates the 
processing of the majer componentS ina digital ccntroller. 
An analog component circuit can be replaced by an analcyg to 
digital converter, a digital processor, anda digital to 
analeg converter. Sore of the benefits which can be realized 
by doing this are; 

it. A high-speed processor could actually prccess a 
number of multiplexed signals, performing processing func- 
tions on a number of independent channels. 

2- The processing function iS permanent in Software, 
unless deliberately changed, and will not drift with age. 

3. The processing function can be changed without 
changing components, merely by changing software. 

4. Accuracy can be made very high and can be changed 
merely by changing scftware. 

5- Processing, which previously required large compo- 
nents such as inductors in low-freguency controllers, can 


now be performed by very small digital circuits. 
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Figure 7.1 ADAPTIVE CONTROLLER 


63 


me 













X(t) DIGITAL Y(t) 
input - PROCESSOR CONVERTER]! output 






er ar ET mS petite 
SS EE 


SS et ee ee eh ee 


i 


Pigure 7.2 DIGITAL BLOCK DIAGRAM 


In converting an analog controller to a digital 
contrcller, the process can be broken down into the 
following steps: . 

1. Determine the desired analog transfer function. 

2- Set the sampling frequency. 

3. Apply the bilinear z-transformation. 

4. Match one point in the s domain to the z domain. 

5e Oktain the optimum constant coefficients. 

6. Obtain the digital transfer function. 

7e Oktain the sisulation diagran. 

The optimal ccntroller parameters can be stored in 
RCROLY. Intel company markets a 4 megabit non-volatile 
read/write bubble semory. It is supported by a VSLI 
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contrcller which provides a black box interface. It is easy 


to use 

SOrS. 
ie 
Ze 
o. 
4. 
De 
6. 
ies 


and can be used with any 8- or 16-bit microproces-~ 
The bubble memcry advantage 1S: 
Fast access time compared with disk or tape. 
Ncen-volatile. 
Wide temperature ranye of operation. 
Werking storace. 
Pcrtable operation 
Low power. 


High reliability. 
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VIII. CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE STU 
Aw CCNCLUSIONS 
In designing the controller, three different ship scdels 
were used. Using the second-order Nomoto model Equation Lak 
allowed comparison of results with Reid's [Ref. 7] [Ref- 10] 
wOrk. It is clear that the answers obtained by function 
Rinimization agree closely with Reid's results as shown in 
Tables 4 and 5. A tketter description of the ship is the 
third-order Nomoto model which involves both the sway and 
yaw eguations. This nrodel includes the two dominating foles 
of the ship. The best model to describe the dynamics of the 
ship is a Taylor's series expansion. This allowes both 
linear ard nonlinear terms in the eguations of moticn to 
affect tke design of the contrcller. | 
Yo determine which controller structure would prcvide 
the minimum cost due to steering, various structures were 
studied. It was found that the dynamics of the plant deter- 
Dines the optimum structure for the controller. In caln 
water study, when using a second-order Nomoto model, the 
best structure was ccntroller A. When the third-order Ncmoto 
model Equation 2.2 was used the best structure was 
contrciler C, but the difference is slight. Observe that in 
each case the contrciler zeros cancel the plant poles. When 
the eguations of motion were used for tne plant, the best 
structure waS contrcller C. When the equations cf notion 
were coupled to a sea state generator and the cost function 
waS fininzized in the presence of a seaway, the best struc- 
ture was ccntroller A. This study concludes that ccntrcller 
A should te used. 
A function minimization subroutine is an engineer's tool 


which can be used in many engineering problems. Previously a 
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Matched filter was designed for the Naval Postgraduate 
Schocl research project on the Space Transportation Systen 
(STS) for the Get Away Special Program. It was natched to 


the signature of the auxillary power unit (APU) on kLoard the 


space shuttle. The goal was to turn on the solid state 
recording system before lift off, to record the acoustic 
power generated inside the shuttle bay. Basically the 


Matched filter is a Finite Impulse Response (FIR) filter 
with the weights calculated to obtain the least squared 
error of the desired output when the input is the signature 
of the AFU. Figure 8.1 shows the scheme used to evaluate 
the FIR weights. 


De SRECORSENDATIONS FCR FUTURE STUDY 

In the future mest ships both military and commercial 
will have GPS receivers as part of their navigation eguip- 
ment. Using extremely accurate three-dimensional fosition 
and velocity infcrmation from satellite platforms will allow 
Ships to navigate accurately in and out of ports. The func- 
tion ginimization Subroutine 1S a powerful Loo 2Og 
desicning the contrcller. This routine simply takes the 
inputs that require finimizaticn and adjusts the parameters 
to acccmyplish this task. The cost function for the added 
drag due to steering is a function of yaw error and rudder 
motion. The use of function minimization and NAVSTAR/GPS 
provides the means for optimization for guidance and cont- 
roll. There are several areas that need future study and 
work. 

1. Should the objective change to track following then 
it 1S necessary to minimize the yaw error only. This would 
be very important both militarily and commercially should a 
port be grined. If the ship could follow a stringent route, 
knowledge of mine locations would allow access. 
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Figure 8.1 MATCHED FILTER DESIGN 


2. A controller for orbit keeping for satellites with 
High-Energy Laser weapons would be very important. The small 
far-field spot size cf a focused laser beam can be selec- 
tively focused on the most vulnerable component on the 
target, Facilitating precision energy deposition, and 
greatly increasing the probability of a kill. 

3. An adaptive controller to minimize track error on 
koard a cruise missile could be programmed for selective 
targets. 
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G- Military and commercial aircraft can benefit just as 


do ships ty reducing drag to minimize fuel consumption. 
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APPENDIX A 
PROGRAM TO CALCULATE OPTIMAL GAILHS 


The proyram iS set up to calculate the optimal gains for 
contrcller A. It 1s referenced in Chapter five and six. Lt 
can easily be modified to obtain optimal gains for the rest 
of the ccntrollers. After obtaining the optimal gains the 
program most be modified to do a Simulation. The prograz has 
sufficient comments for appropriate changes. It is refer- 
enced in Chapter two. 

This program can be modified to obtain the Ncneto 
models. It is referenced in Chapter two. The following need 
to be changed. 

C GAIN COSFFICIENTS [0 Bree ere oD 
K=XX (1) : 
TP T=KK (2) 
TZ=XX (3) 
TP2=XX (4) 
ERRCK SIGNAL TO LEIVE RUDDER (YAW ACTUAL —- YAW COMMANT) 
FCR ECUATIONS OF FOTION. 
D=YAW ~- YAWC | 
C ERROR SIGNAL TO CRIVE RUDDER (YAW COMMAND —- YAW ACTUAT) 
C FCR NCMOTO 35D ORIER MODEL. 
D2=YAWC-YAW2 
X1= (D2-X2)/TP1 
X¥3=K* (TZ*X1+ XZ) 
X4= (X3-X5) /TE2 
C INTEGRATION 
X2=X2+X1*DELT 
X5=X5+X4*DELT 
YAW2Z=YAW2+X5*DELT 
C CCST FUNCTION 
IDIFF=TDIFF+ (YAW-YAW2) **2 
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2' FT YDOT=",F8.4,' FI/SEC',/,2x,' UC=',F8.4, 
3' FT/SEC U=',F8-4,' FIVSEC ° UboT=',F10.6, 
HoY FT/SEC**2 v=',Fs8.4,' FT/SEC vVDOT=",F10.6, 
5! FT/SEC#*2" ,/,2X,4YAWC=",FS.4," DEG YAW=",F15.7, 
6' DEG YAW RATE=',715.7,' DEG/SEC YAW ACCEL=! 
7,F15.7,' DEG/SEC*¥*21,/,2X," RUDDER =',F15.7," DEG ',/) 
TCCUNT=1 
C TEST IF WANT TO STOP 
300 IF (TINE.GE.ETINE) GO TO 400 
C INTEGRATION STEP SIZE DELT 
DELT=1.0 
C IKIEGRATION 
U=U+UDOT*DELT 
V=V+VDOT*DELT 
R=R+RDOT*DELT 
YAW=YAW+R¥*DELT 
X 2=X2+DK2*DELT 
C CCNVERT SHIP TO FIXED COORDINATES ON EARTH 
EDOT-U*DCOS (YAR) “VADSIN (28) 
YDOT=U*DSIN (YAK) +V*DCOS (YAW 
X=X+XDOT* DELT 
Y=Y+YDOT*DELT 
TIME={TIME+DELT 
ICCUNT=ICOUNT#1 
ISE=ISE + LAMCA*YAWE**2 
ISE=ISR + D*®*2 
GG TO 200 
J= TDIFF = COST FUNCTICN 
40C  ‘IDIFF=ISE+ISR 
WRITE (6 (500) ISE,ISR,TDIFF,K1,T1,T2 
500 FORMAT ( i ip bowel F15.7,' ISR=',F15.7,' TOTAL=! 
1,F715.7,2X,¢K1=",F 1507, 2%, 'T1=",F 15. 7,2X,'T2=", F156 
REWIND 12 
RETURN 
EX 
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DELETE ALL THE FOCILOWING SUBROUTINE IF SIMULATIONS Ghay, 
AND NOT OPTIMIZATICN iS WANTED 
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SUBROUTINE BOXELX (CATEGORY HO) 

PURE Oss 
BCXPLX IS A SUBROUTINE USED TO SOLVE THE PROBLEMA CF 
leocacting A MINIMUF (OR MAXIMUM) OF AN ARBITRARY OBJECT- 
ive function SUBJECT TO ARBITRARY EXPLICIT AND/CR 
implicit constraints by tHE COMPLEX METHOD Ot =h QU See 
explicit constraints are dEFINED AS UPPER AND LOWER | 
bounds on the independent variables IMPLICIT constraints 
May be arbitrary function of the varIABLES. TWO FUN- 


ction subprogram tc evaluate the objective FUNCTION AND 
ee ae con PERT tg otc” coe must be SUPPLIED 


y the user (see EXAMPLE BELOW). BOXPLX ALSO HAS tHE 
Shae to perform integer DEG ena 0 Pee the values 
of the indeperdent variables are restricted to integers. 


USAGE 
CALL BOXPLX (NV,NAV,NPR,NTA,R,XS,1IP,XU,XL,YMN,IER) 


DESCRIPTION OF PARAMETERS 


NV AN INTEGER INEUT DEFINING THE NUMBER OF INDEEENDT 
VARIABLES OF THE OBJECTIVE FUNCTION TO BE MININI 
NOTE: MAXIMUM NV + NAV IS PRESENTLY 50. MAXIMIM N 
25. IF THESE LIMITS MUST BE EXCEEDED, PUNCH A SO 
DECK IN THE USUAL MANNER, AND CHANGE THE DIMEN 
STATEMENTS.~ 
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NAV AN INTEGER INEUT DEFINING THE NUMBER OCF AUXILIARY var 

LaELES THE USER WISHES TO DEFINE FOR HIS OWN CONVENIENCE. 

TYPICALLY HE ABY MSH TO DFFINE PHE YALUE OF “BAGH IMELICI 

CONSTRAINT FUNCTION AS AN AUXILIARY VARIABLE. beet HLS 

IS DGCNE, THE OPTICNAL OUTPUT FEATURE OF BOXPLX CAN BE 

USED TO OBSERVE THE VALUES OF THOSE CONSTRAINTS AS THE 

SCLUTICN PROGRESSES. AUXILIARY VARIABLES, IF USED, 

SHCULD BE EVALUATED IN FUNCTION KE (DEFINED BELOW). 

NAV MAY BE ZERQ.,. 

NPR INFUT INTEGER CONTROLLING THE FREQUENCY OF OUTPUT 

aqgeseired for ede nesetic Poe cses. 

IF NPE .LE. QO QO CUTPTU EL. SE 

PeOoOlcuceD BY BOXPLY - OTHERWISE, THE CURRENT COMPLEX CF 

K= 2*KNV VERTICES ARD fe eer CENTROID WILL BE OUTPUT AFTER 
EACH NPR PERMISSIBLE TRIALS. Tok eVUMBER OF TORAL TREATS, 
pombe OF FRASTSBLE TREALS, NUMBER OF FUNCTION EVALUATICNS 

AND NUMBER OF IMPLICIT CONSTRAINT EVALUATIONS ARE IN- 

CLUDED IN THE OUTEUT. 

ADDITICNALLY, (WHEN NPR .GT 0) THE SAME INFORMATION 

WELL BE CUTPUT: 

ipeeheetHe DNITEAL FOINT IS NOT FEASTIELE, 

2 RRFIER THE FIRSEZ COBPLETE COMPLEX IS GENERATED, 

3 RPOAP Er raASLELE VERTEA CANNOT BE FOUND AT SOME TRIAL, 

4 IF THE OBJECTIVE VALUE OF A VERTEX CANNOT BE MALE 

NO~LCNGEER-WORSTI. 

2} Treas Lee itT ON TRIALS NTA) IS REACHED AND, 

6) WHEN THE OBJECTIVE FUNCTION HAS BEEN UNCHANGED FOR 

2*NV TRIALS, INDICATING A LCCAL MINIMUM HAS BEEN 

FOUND. 

Mmreteer USER WHESHES TO TRACE THE PROGRESS OF A SOLUTICN, 

A CHOICE OF NPR = 25, 50 OR 100 IS RECOMMENDED. 

NTA INTEGER INPUT OF LIMIT ON THE NUMBER OF TRIALS 

allowed in the calculation. 

PeelHE USER INPUTS NTA .LE.- A default 

VALUE CF 2000 IS USED. WHEN "Suts LIMIT IS REACHED 

CONTECL RETURNS TO THE CALLING PROGRAM WITH THE BEST 

ATTAINED OBJECTIVE FUNCTION VALUE IN YUN, AND THE BEST 

ATTAINED SOLUTION ECINT IN XS. 

R A REAL NUMBER INPUT TO DEFINE THE FIRST RANDCM NGMEER 
USED IN DEVELOPING THE INITIAL COMPLEX OF 2¥*NV VERTICIES. 

ne tes RG iota, . 1% IF IS NOT @iTHIN .THESE EOUNCS, 
Beeeen BE REPLACEIL BY 1./3. -« 

KS INPUT REAL ARHAY DIMENSIONED AT LEAST NV+tNAV. 

the first nv must contain a 

FEASIBLE ORIGIN FCE STARTING THE CAL- 

CULATICN. THE LAST NAV NEED NOT BE INITIALIZED. DECN 

RETUEN PROM BOXPLX, THE FIRST NV ELEMENTS OF THE ARRAY 

CONTAIN THE COORDINATES OF THE MINIMUM OBJECTIVE 

function, AND THE REMAINING NAV (NAV .GE. QO CONTAIN THe 

values or THE CORRESPONDING AUXILIARY VARIABLES. 

IP INTEGER INPUT FOR OPTICNAL INTEGER PROGRAMMING. 

if 1lp=1, THE VALUES OF THE INDEPENDENT VARIABLES WILL 

be replace ed WITH INTEGER VALUES (STILL STCRED AS REAL*4). 

XU A REAL ARRAY CIMENSICNED AT LEAST NV INPUTTING THE 

upper BCUND ON EACE INDEPENDENT VARIABLE See EXPLICIT 

CONSTRAINT). INPUT VALUES ARE SLIGHTLY ALTRRED BY BOXPLX. 

XL A REAL ARRAY CIMENSIONED AT LEAST NV INPUTTING THE 

lower bound on each independent 

VARIABLE, (EACH EXELICIT CONStraint). 
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NOTE; FCR BOTH XU AND XL CHOOSE REASONABLE 
VALUES IF NONE ARE GIVEN, NCT VALUES WHICH ARE 
Dagnitudes ABOVE CR BELOW THE EXPECTED SOLUTION. 
1nput values are SLIGHTLY ALTERED BY BOXPLX. 
YMN THIS OUTPUT IS THE VALUE (REAL*4) OF THE OBJECTIVE 
funcTICN,CORRESPONLING TO PHE SOLUTION POENT OUTPUT “Tikes 
IER INTEGER ERROXH RETURN. TO BE INTERROGATED UECN 
return FRCM BCXPLX. IER WILL BE ONE OF THE FOLLOWING: 
=~} CANNOT FIND FEASIBLE VERTEX OR FEASIBLE CENTROID 
AT THE Sine OF 4 RESTA KSEE ‘METHOD’ BELOW). 
= FUNCTION VALUE UNCHANGED FOR ‘N* TRIALS. (WHERE 
THIS IS THE NORMAL RETURN PARAMETER. 
=1 CANNOT DEVELCE FEASIBLE VERTEX. 
=2 CANNOT DEVELOEF A NO-LCNGER-WORST VERTEX. 
=3 LIMIT ON TRIALS REACHED. NTA Sean, 
NCTE: VALID RESULTS MAY BE RETURNED IN ANY OF THE 
ABOVE CASES. 

EXAMELE CF USAGE 
THIS EXAMPLE MINISIZES THE OBJECTIVE FUNCTION SHOWN IN 
the EXTEENAL FUNCTICN Se ee THERE ARE TWO INDEPENDENT 
Var tatres ae & X(2), AND TWO IMPLICIT CONSTRAINT 
funcizon xf ) & X{4) WHICH ARE EVALUATED AS AUXILIARY 
variatles (see EXTERNAL FUNCTION KE({X) ). 


DIMENSICN XS(4),XU(2),XL(2) 
STARTING GUESS 
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FUNCTION KE ee 
EVALUATE CONSTRAINIS. SET KE=0 IF NO IMPLICIT CONSTRAINT 
1s vViCLATED, OR SET KE=1 IF ANY IMPLICIT 
ccnstraint 1S violated. 
DIMENSION X (4) 


X1 = X 

X2 = x2 

KE = 0 

X(3) = X¥1 4+ 1.732051*x2 

TE (X(3) LT, 9. OR. X(3) .GT. 6.) GO TO 1 
X (8) = 81/1. 7352051 =x 
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TEth) .GEs C.) /RETURN 
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METHOD 


EOGOMPLEX METHOL IS AN EXTENSION AND ADAPTION OF 
€ simple method cf linear le poiat an 
ARTING WITH ANY CNE feasible point in n-dimension 
NCCMPLEX" OF 2*N vertices 1S constructed by 
LECTING RANDOM ECINTS WITHIN THE feasible 
GICN. FOR THIS FURPOSE N COORDINATES ARE FIRST 
NDC¥LY CHOSEN WITHIN THE SPACE BOUNDED BY EXPLICIT CCN- 
Peo Pomel DEFINES A TRIAL INITIAL VERTEX. 
s then checked for possible violation 
MEFLICIT CONSTRAINTS. IF one or more are violated, 
HIAL ENITIAL VERTEX IS DISPLACED half of its 
Ch Pues THteCeENieOlp OF PREVIOUSLY SELECTED initial 
Poe eee be tosAns THIS PISPLACEMENT PROGESS IS 
2 UNG THE VE hao DSCOMbessAStbLs. IF THIS 
TO ageee afte +10 displacements, 
N AFTER EACH VERTEX IS ADDED 

centroid is checked fcr 
ASLELE, gene -last trail 
EFFORT TO GENERATE an alter- 

I A IF 5*N+10 VERTICES ARE 
DCNED CONSECUTIVELY, THE SOLUTION IS TERMINATED. 


N INITIAL COMPLEX IS ESTABLISHED, THE BASIC 

utation loop is initiated. 

E INSTRUCTIONS FIND THE CURRENT WORST vertex, that 
HE VERTEX WITH THE LARGEST CORRESPONDING value for 
BJECTIVE FUNCTION, AND REPLACE THAT VERTEX BY 
IVER-REFLECTICN THROUGH THe comolaoLrpeor mtlL CHHER 
ces. 1f the vertex to be 

ae IS CONSIDERED AS A VECTO -~space, 
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ER-REFLECTICN IS OPPOSITE TH DIRECTION, IN- 

D IN LENGTH EY THE FACTOR 1.3, AND COLLINEAR WITH 

PLACED VERTEX AND CENTROID OF ALL OTHER VERTICES.) 
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AN OVER~REFLECTION IS NOT FE 
a7 IS CONSIBERED NOT-PERMNIS 
ISPLACED EALFWAY TOWARD THE 


1 IS MADE BY REFLECT 
RESENT BEST 

ROUGH THE CENTROID. IF 5¥*n+10 
LECTIONS OCCUR WITHOUT A 
ESULT, THE CURRENT BEST 
i E 
E 
K 


Ur 


wna 


SSeS I CIO) eS 


wo b+ 
eS tie = OO Sito bts 


aoe POINT FOR A 


ES 

EN "S env+10 CONSECUTIVE 
GNIFICANT CHANGE IN THE 
A 
T 


“Mrgy we wmNOm 
QM wWNORUUNHID > 
= Ord sitttd mre 


bet bed Oo ta oe 


: 
. «tN ALL CASES 

LAST RESTART DID NOT 
NT IN THE MINIMUM 


HQHMAG’a thy abory 
W 
rye tthe Oty 
Ht We Ho Om cs 
OMHOROnmiNIN Mrs 
toe Fr + fe 
<—Qheqrazri ry 
AammHOmmot: 


A 
h 
iL 
E 
E 
A 
aL 
Lt 
FE 
A 
if 
E 
E 
A. 
HE 
Cc 
E 
E 
my 
B 
i 
oO 
Hi 
5 
HE 
a 
S 
E 
I 
R 
H 
H 
0 
ND 
a 
V 
i 
E 
I 
U 
E 
ze 
E 
K 
A 
E 
K 
- 


mgt etnodedine Pate HPOWHOURMHHHNAQH Peder SC OHOP bbe NicthHs 


be ed ee Db bn He ed i eH 
= Om Woops WL) be bd bd 
MMO WH YOs Athy 

O Mtoe Heth WHO 
ore mm lctsHrdid 


wl 


ANNA AQAAA AANA|A QNQNAQAANAAA AQANQQAQQAAQAAQAQAAAAA QQAQAQA AAA AQAAAA AQAA AAANAAQAAAQAAAANA 


IT IS RECOMMENDED THAT THE USER READ THE REFERENCE FOR 
FURTHER USEFUL INFORMATION. IT SHOULD BE NOTED THAT THE 
ALGORITHM DEFINED THERE HAS BEEN ALTERED TO FIND THE 
CONSTRAINED MINIMUMS, RATHER THAN THE MAXIMUM. 
REMARKS 

THE INTEGER PROGRAMMING OPTION WAS ADDED TO THIS EROGRAM 
AS SUGGESTED IN REFERENCE (2). BA Seep 
TRE SL ACO ae variable version of DOX lx 
WOULD BE EASY TO CREATE BY DEClaring "1p" to Dewag arcay 
OF NV CCNTROL VARIABLES WHERE IP (1)=1 would indicate 
THAT THE I-TH VARIABLE IS TC BE CONFINED to integer 
VALUES. BACH STATEMENT OF THE FORM “LF (iP . EGC che 
WOULD THEN NEED TC BE ALTERED TO "ZF (1P(1).. EO. 1)" Giese 

WHERE THE SUBSCRIPT IS APPROPRIATELY CHOSEN. NCRMALLY, 
fu AND XL VASES ARE ALTERED TO BE AN EPSILON "WITHIN! 
actual valu 
DECLARED BY. THE USER. THIS ADJUSTMENT IS NOT MADE 
WHEN IP=1. 
NCTE: NO NON-LINEAR PROGRAMMING ALGORITHM CAN GUARANTEE 
that the answer found is the global 
MINIMUM, RATHER THAN JUST A local Minimum. however, 
ACCORDING TO REF.2Z, THE COMELEX method has an advantage 
-IN THAT IT TENDS I0 FIND THE GLOBAL nDinimum more 
FREQUENTLY THAN MANY OTHER NCN-LINEAR PROGRAM - 
MING ALGORITHMS. 
IT SHCULD BE NOTEL THAT THE AUXILIARY VARIABLE FEATURE 
CAN ALSO BE USED IC DEAL SWiT8 
PRCELEMS CONTAINING EQUALITY CONStraints. any eguality 
CCNSTRAINT IMPLIES THAT A GIVEN VARiable is not truly 
INDEPENDENT. THEREFORE, IN GENERAL, ONE variable 
INVCLVED IN AN EQUALITY CONSTRAINT CAN BE RENUMBERED from 
THE SET -OF NV INDEEENDENT VARIABLES AND ADDED TO THE SET 
OF NAV AUXILIARY VARIABLES. THIS USUALLY INVCLVES 


Fe Due THE INDEPENDENT VARIABLES OF THE GIVEN 
Le 
gUBROUTINES AND FUSCTICNS REQUIRED 


SUBROUTINE 'BOUT' AND FUNCTION *FBV' ARE INTEGRAL 
parts of THE BOXPLX PACKAGE. 


THO FUNCTIONS BUST BE SUPPLIED BY THE USER: THE FIRSIe 
ne A) is used to evaluate the implicit 

C RAINTS. SET KE=0 AT THE beginning of the function 
N EVALUATE THE IMPLICIT CONStraints. in the exazple 
Pee TG eee ae CON eee pear hust be within the 

« « LEcadk (3) See oe E SECOND constraint x (4) 
E .~GE. QO. IF ELT ER CONSTRAINT IS not wit aoe 
UNDS CONTEOL IS TRANSFERRED TO STATEMENT 

S SET TO "1" AND CONTROL IS RETURNED TO BOXSIx. 


COND FUNCTICN se USE MUST PROVIDE EVALUATES THE 
jective function. it 
LED eA AS SHCWN IN THE Ee above, and fe 

os T TC THE VALUE OF THE OBJECTIVE function 
RESPCNDING TO CURRENT VALUES OF THE NV INDEPENDENT 
IAELES IN ARRAY 'X*, 


REFERENCES 
BOX, M. Je, "A NEW METHOD OF CONSTRAINED CPIiMi Zari 
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SUBRCUTINE BOXPLX (NV,NAV,NPR,NT2Z,RZ,45,1P,8B0,BL, (MN ,LTER) 


DIMENSION V (50,50 FUN (50 SUM (25 CENGZ> XS(NV 
viuinyy DLINY)” ys (20), (25), (25), (NV) 


Ay = 95 

EEF = 1.E-6 

NIA = 2000 

IF _(NTZ.GT.0) MIA = NTZ 

kK = 

Pm (ieewoeon. ORSROGE. 1. )R=1.73. 
NVT = NV+NAV 


Honebew ARO, EXPLICIT PLUS IMPLICIT 

” CURRENT TRIAL NO. 

z ~ CURRENT NO. OF PERMISSIBLE TRIALS 
Mae ORRENT NO. OF TIMES F HAS BEEN ALMOST UNCHANGETL 
CHECK rEASIBWeCL TY OF START POINT 


=z 
i 


DC 4 I=1,NV 
VI = XS ( ) 
IF (BL { ).LE.¥1) GO TO 1 
ie = = 
VI = BL (I) 
GC TO 2 
1 IF ooo GO TO 3 
Tl 
VI = Bh 
2 IF (NDR. G 0) KRITE (6,49) II 
3 V(t, = 
Cin) = vt 
IF {1P.£Q.1) GC TO 4. 
BL =P Br hy eam axl EP *ABS (BL (I 
BU{I) = BU(T)-AMAX1 EP*ABS (BU (I 
4 SUM(I) = V 
NCE = 
NUMBER CF CONSTRAINT EVALUATIONS 
V{1,1))-EQ.0) GO TO 5 
IF { vite! {6 TO. 12 
WRITE {6 ,50) 
GOuie 12 
S NFE = 1 


KUREER OF VERTICES (K) = 2 TIMES NC. OF VARIABLES. 
K = 2*NY 
NUMBER OF DISPLACEMENTS ALLOWED. 
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NLIM = S¥*NV+10 


NUSBER OF CONSECUTIVE TRIALS WITH UNCHANGED FE TC 
tecrninate. 

NCT = NLIM+NV 

ALPHA = 1.3 


FK = K 
FKM = FK-1. 
BETA = ALPHA+1. 
INSURE _SZED OF RANDOM NUMBER GENERATOR IS ODD. 
IF (MOD {ZOR,2)-EQ.0) ICF=IQR+101 
SET UP INITIAL VERTICES 
FUN{1) = FE(V (1, 1) ) 
YUN = FUN (1) 
6 FI = 1. 
FUNCLD = FUN(1) 
DC 15 I=2,K 
FI = FItil. 
LIMT = 0 
7 LIMT = LIMT+1 
END CALCULATION IF FEASIBLE CENTROID CANNOT BE FOUND. 
IF (LIMT.GE.NLIIM) GO TO 11 


DC 8 J=1,NV 
BANCDCM NUMBER GENERATOR (RANDU) 
ICR = IQR *6553S 


CG 
IF (IQR=LT.0) IQR = IQR+2147483647+1 


13E- 
V(J,I) = BL(J)+ROX*(BU(J)-BL(J 
(Shs oT) SALT (0 fT). 5) 


La EOwOe GCONTemS 


Mock 


(rb. E 


} VI= AINT(VT+.5) 
S Vidgl)"= 


v 
(J, ,1)+CEN(J)) 
T 


10 CCNTINUE 
11 IF (NPR. LE. 0) [Go Bomar 


TE 
Li adyt (it. NPT, NFE,NCE,NV,NVT,V,1,FUN,CEN,1) 
48 
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ASIBLE CENTROID FOR STARTING. 


E 
C.0) GO TO 60 
J) -V¥(J,1) 


tot = tr} 
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END CF LOOP eee OF INGZTIAL COMPLEX. 
IF (NPR.LE. Wh, RE HOR U7 
CALL BOUT Dik Be oN ,nivi,V¥,K,rUN,CEN, 0) 


FIND ae WORST VERTEX, TEE ‘'J'TH. 


DC 16 I=2,K 
TE (EON (JJ. GE. FUN (1) ) Gc Tome 


16 CCNTINUE 


Cy 


BASIC LOOP. ELIMINATE EACH WORST VERTEX IN TURN. 
it twust become NC LONGER WORST, NOT MERELY IMPROVED. 
find next-to-vertex, THE 'IN'TH ONE. 


17 JN = 
faenonty dN = 2 


fre ~EO.J 

FUN (JN). 

18 CONTINUE 
LIMIT = NUMBER CF MOVES DUBING THIS TRIAL TOWARD THE 
centroid, DUE TO PUNCTION VALUE. 


CCHEUTE CENTROII AND OVER REFLECT WORST VERTEX. 


GG Tom 
GE.FUN(I)}) GO TO 18 


VI= V J) 

Rays ="SUM(I)-VT 

CEN (I) = SUM{I)/FKM 

VI = BETA*CEN {I) -ALPHA*VT 

Pm (1P.fO<1) VI = NENT (VT+.5) 
INSURE THE EXPLICIT CONSTRAINTS ARE OBSERVED. 
19 V(I,J) = AMAX1(AMIN1(VT,BU (TI) ),BL(T)) 

NI = NT+1 


CHECK FOR IMPLICIT CONSTRAINT VIOLATION. 
20 De 25 N=1 ¢ NLIM 
NCE +1 
IF (RE(Y (1, By) .50. 0) Gor ro 26 


EVERY 'KV'TH TIME, OVER-REFLECT THE OFFENDING VERTEX 
thrcouch the BEST (ERTEX. 


IF fRgD (me KV PRE.O0) GO TO 22 
CALL K,FUN,M) 
DC 21 I=1,NV 
VI = BETA*V (1,®) -ALPHA*V (I,J 
MEMGEE. E021) V1 = AINE (V +25 
21 V{I,dJ) = AMAX1(AMIN1T(VT,BU(L)) ,BL(L)) 
GC TO 24 
CCNSTRAINT VIOLATION: MCVE NEW POINT TOWARD CENTROID. 
22 DC 23 I=1 
VI = SEE IEG J) 
IF (IP v ALN (VT+. 5) 
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24 NT = NT+1 
25 CCNTINUE 
IER = 1 
CANNOT GET FEASIELE VERTEX BY MOVING TOWARD CENTROID, 
CR EY OVER-REFLECTING THRU THE BEST VERTEX. 
IF (NPR.LE.0) GO TO 42 
WRITE (6,5 ) NI nd 
CALL BO BOUT (NT,NDT,NFE,NCE,NV,NVT,V,K,FUN,CEN,J) 
G 
FEASIBLE VERTEX FOUND, EVALUATE THE OBJECTIVE FUNCTICN. 
€ NEE = NFE# 
FUNTRY = PE(V({1,d)) . 
TEST TO SEE IF FUNCTION VALUE HAS NOT CHANGED. 
AFQ = ABS (FUNTRY-FUNOLD) 
AMX = ANBYT (AES (ED¥FU NO D) , EP) 


ACTIVATE THE FOLIOWING TWC STATEMENTS FOR DIAGNCSTIC 
Furrposes onsg; | 
A Rae we 99) J,AFO,ANX,FUNTRY,FUNOLD, FUN (J), FUN (JN) 


99 “EORAR (1X,13,6E15.7, 215) 
E (APO. GT: AMX) GO TO 27 


NTE TFS+1 
Iz (NTPS-LT.NCI) GO TC 28 
IF (NPR2LE.O) GO TO 42 
WRITE é3 3) 8 
CALL B Bh (NT ,NPT, NFE,NCE,NV, NVZ,V,XK,FUN,CEN, 0) 
GO TO 42 
27 NIFS = 0 
IS IEE NEW VERTEX NC LCNGER WORST? 
28 IF (FUNTRY.LT.FUN(JN)) GO TO 34 
TRIAL VERTEX IS STILL WORST; ADJUST TOWARD CENTROID. 
EVERY 'XV'TH TIME, OVER-REFLECT THE OFFENDING VERTEX 
through the BEST VERTEX. 
LIM? = LIMT+1 
IF (MOD (LIMT, KV) -NE.0) GO 70 30 
CALL FBV (K,fUX, i) 
DC 29 T=1,NY 
VI = BETA*V (I Ht) -ALPHA*V (I,J 
IF (IP,EQ.1) 1 = AINT (V +25 
29 V(I,J) = AMAX1{(AMIN1(VT,BU (Z)),BL (I)) 
GC TO 32 
30 DO 31 


VI = 5+ \eEN I)+V 
IF (IP. 20. 1) AEG ALNTUVT+.5) 
Ve: o yr! 

31 CONTINUE 

32 IF (LIMT.LT.NLIM) GO TO 33 


CANNOT MAKE THE tJ'TH VERTEX NO LONGER WORST BEY 
dGisrlaciny toward 


THE = a OR BY OVER-REFLECTING THRU THE BEST VERTEX. 


TF NPR -LE- 0) GQ To 42 
AL “Boge” <k 
ALL | ( 1,xPf, NFE, NCE,NV,NVI,V,K,FUN,CEN,J) 
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GC TO 20 
AMROPLACKMENT FOR VERZE J. 


NET = NPT+ 
FVERY 100'TH PERSIS 
summation to AVCIL 

IF (HOD (NPT,100C). 


DO 36 I=1,NV 
Sumi) = 0- 


DC 35 N=1,K 
M(I) = SUM(I)+V(I,QN) 


M (I 
N(Z) = SUM (I) /PK 


LE TRIAL, RECOMPUTE CENTROID 
EPING ERROR. 
POnconi0 37 
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39 IF NBR«LE=0) GO TO 
IF (MOD (NPT, PE) «NE. 0) GOmronnG 


CALL BOUT (NT, ,NPI,NFE,NCE,NV,NVI,V,K,FUN,CEN, IC) 


FAS THE HAX. NUBEER OF TRIALS BEEN REACHED WITHCOT 
pew erence ¢ Tresor, GO TO Now TRIAL. 
40 IF (NT.~GE.NTA) Co. TO 41 


ee CeO RST VERTEX NOW BECOMES WORST. 


GC TO 17 
41 IER = 3 
IF (NPR.GI.9) HRITE (6,54) 


SUrerecrOR POINT FOR ALL ENDIN 
CANNOT DEVELCE FEASIBLE VE 
CANNOT DEVE A NO-LCNGER- 
FUNCTION VA 
bee PON Th 
CANNOT FIND 

42Z CCNTINUE 
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TF (IER. ‘Gh. 3) "66 ko 44 
RESTART IF THIS SOLUTION IS SIGNIFICANTLY BETTER THAN 
cue PERL OUSs OR lf eo tomers rio. TRY. 
NPR LE: 29) GO TO 43 
WEI (6,5 (U,YON, PON (1) 
43 Fu (A): E- YAN) GO 
TF (ABS {F N(M)-YHN). LE. AMAX1 (EP, EP* YM) ) Gor Powe 7 
GIVE IT ANOTHER TRY UNLESS LIMIT ON TRIALS REACHED. 
44 oYMN = FUN (4) 
FON(1) = FUN (M) 
DO 45 I=1,NV 
CEN s = x (zoe 
Sum(I) = vV{I,M 
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u5 Vil,1) = V(z,%) 
DC 46 I=1,NVT 

4ué XS(I) = V(I,M) 
IE ea hae GO TO 6 

4u7 IF (NPR.LE.O) GO TO 48 
CALL BOUT An ,NPT,NFE,NCE,NV,NVI,V,K,FUN,V(1,M) ,-1) 
WRITE (6,56) FUN (M) 

Y& RETURN 


49 FORMAT (200° DE st? DIRECTION OF OUTLYING 
Tvariable at start 
50 FORIAT sO ho TsBLict CONSTRAINT VIOLATED AT 
star €a end 
S51 FCEMAT ('OCANNCT IND FEASIBLE',14,'TH VERTEX SOR 
IceEntroid at. start.:' 
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TI 
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bay DO 
roe 
= tO 
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Ut 
Nn 
Peace hte ley 
xj 


De 1 I=2 
TF _ (ro (4). -LE.FUN(I)) GO TO 1 


1 CCNTINUE 


RETURN 
END 
SUBROCUT 
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same 


NE BOUT ‘ms 
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if 
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Omari bh 
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NV 
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3 E 
(EK ONE 0) =GOn Lom 
(6,7) oO I=1,NV) 


N 
Ke GE. 39) 3 
(6,8) (Cit), OIE 1 ,NV) 


B 
WRITE (6,9) IK, (C(I) ,1=1,NV 
Sa r) e{C (1) , T= 1,NV) 


4 FCRMAT (*ONO. TOTAL TRIALS = ',15,4X, 


I oo MOE 


A) 

re 

biteitt = he 
C3. CH 


oY) 


1'no. feasible trails = 15, 4x 

2'NO. FUNCTION EVALUATIONS =e oe 

3'nc. constraint evaluations = ',1 

ye FUNCTICN VALUE',6X, "INDEPENDENT VARIABLES/ 
SdependENT OR IMPLICIT CONSTRAINTS! 
S FCRMAT (1H ,E18.7,2X, 7E1.7/(21X, 7214 ~ 

6 FORMAT {21X,7E10.7) 

7 FORMAT (1OHOCENTROZD 11X,7E14.7/(21X g7Et ™)) 

8 FCRMAT ('0 BEST VERTEX',.7X,7E14.7/(21X, 1B 7)) 


EN 
Peo tolh UD e 


Z 
J//GC.FTI2ZFO001 DD DISF=SHR, DSN=MSS.52160.A341 
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APPENDIX B 
EXABSPLE PROBLEM USING ICSOS 


The purpose of this example is to demonstrate the 
performance of the prcgram. Consider the control systen of 
Figure 4.1 with controller C. Figure B.1 shows the block 


diagram to evaluate the contrcller parameters. 


-_* 
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The differential equations describiny the system and 


desired performance are: 
¥2=DX2 : 
X4=DX4 
XS=DX5 
R =DR 
YAW=R 

Defining the following cost function: 
a= f (LanpasyawE+42+D*#2) at 

Defining the special functions: 
YAWE=YAWC-YAW 
DX2= (YAWE-X2) /12 
X3=K1* (T1¥*DX2+X2) 
DX4= (X3-X4) /T4 
d= (t3¥dx4+x4) 
dx5=(d-x5)/tp1 
x@=k* (tz *dx5+4x5) 
dr= (x8-r) stp2 

Defining the constants: 

YAWC=1.0 
K=. 14771 
T2=11.2833 
TF1=6.4699 
TE2—5o57 931 
LAMDA=4.2 


and using YAWC=1.0 the optimal solution found by 


the frograr is: 
K 1=.4179916 
T 1=53.69932 
T2=4.970023 
T3I=6.294369 
T4=13.85919 
COST=68. 04735 
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its 


Table 34 shows the specifications of this problem with 
the free parameter cptimum values found. Figure B.2 shows 
the actual yaw and rudder response. 
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